home *** CD-ROM | disk | FTP | other *** search
/ InfoMagic Internet Tools 1993 July / Internet Tools.iso / RockRidge / info-service / prospero / PRM / src / testprog / ocean1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-04-07  |  23.0 KB  |  710 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. #define MAIN_PROG
  5. #include "comm.h"
  6. #include "ocean.h"
  7. #include "CMMDlib.h"
  8.  
  9. char *progname;
  10.  
  11. main(argc, argv)
  12. int argc;
  13. char **argv;
  14.  
  15.    {void Ocean_Work (), Ocean_Idle (), Compute_My_Coordinates ();
  16.     int num_iter;
  17.  
  18.     init_task(argv);
  19.  
  20.    CMMD_receive_bc_from_host ((char *) (& num_iter), sizeof (int));
  21.     fprintf(stderr, "(task %d) received bcast. num_iter = %d\n", _my_taskid, 
  22.         num_iter);
  23.     Compute_My_Coordinates ();
  24.  
  25.     if (my_node < TOT_PROCS)
  26.        Ocean_Work (num_iter);
  27.     else
  28.        Ocean_Idle (num_iter);
  29.     exit(0);
  30. /*    CMMD_sync_with_host (); */
  31.   }
  32.  
  33.  
  34. void Ocean_Idle (num_iter)
  35. int num_iter;
  36.  
  37.    {int inx, all_flag;
  38.  
  39.     do {all_flag = CMMD_reduce_int (0, CMMD_combiner_ior); 
  40. /*        if (my_node == 100) printf ("one: %d\n", all_flag); */
  41.        } while (all_flag == TRUE);
  42.  
  43.     (void) CMMD_reduce_double (0.0, CMMD_combiner_add);
  44.  
  45.     for (inx = 0;  inx < num_iter;  inx ++)
  46.        {do {all_flag = CMMD_reduce_int (0, CMMD_combiner_ior); 
  47. /*            if (my_node == 100) printf ("two: %d\n", all_flag); */
  48.            } while (all_flag);
  49.  
  50.         (void) CMMD_reduce_double (0.0, CMMD_combiner_add);
  51.  
  52.         do {all_flag = CMMD_reduce_int (0, CMMD_combiner_ior);
  53. /*            if (my_node == 100) printf ("three: %d\n", all_flag); */
  54.            } while (all_flag);
  55.        }
  56.    }
  57.  
  58.  
  59. void Ocean_Work (num_iter)
  60. int num_iter;
  61.  
  62.    {void Compute_Constants (), Initial_Matrix (), Exchange ();
  63.     void Compute_Tau (), Do_One_Iteration (), Compute_F ();
  64.     void Solve_Helmholtz (), Print_Mat ();
  65.     int inx, jnx, knx;
  66.     double partial_sum, time, Compute_Sum ();
  67.  
  68.     only_on (OBSERVER) Second_Set (N_TIMER);
  69.     only_on (OBSERVER) 
  70.         {time = Second_Look (N_TIMER);
  71.          CMMD_send (CMMD_host_node, TIME_TAG, &time, sizeof (double));
  72.         }
  73.  
  74.     Compute_Constants ();
  75.  
  76.     Compute_F ();
  77.  
  78.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  79.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  80.             Psi_Um[inx][jnx] = 0.0;
  81. /*    only_on (OBSERVERS) Print_Mat ("Psi_Um", Psi_Um);  */
  82.  
  83.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  84.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  85.             Psi_Lm[inx][jnx] = 0.0;
  86. /*    only_on (OBSERVERS) Print_Mat ("Psi_Lm", Psi_Lm);  */
  87.  
  88.     Initial_Matrix (Psi_B,  1.0, 1.0, 0.0);
  89.     Exchange (Psi_B);
  90. /*    only_on (OBSERVERS) Print_Mat ("Psi_B", Psi_B);  */
  91.   
  92.     Solve_Helmholtz (eig_2, Psi_B, Psi_B);
  93. /*    only_on (OBSERVERS) Print_Mat ("Psi_B", Psi_B);  */
  94.  
  95.     partial_sum = Compute_Sum (Psi_B, 0.25, 0.50, 1.00);
  96.     psi_bi = CMMD_reduce_double (partial_sum, CMMD_combiner_add);
  97. /*    only_on (OBSERVERS) printf ("psi_bi = %e\n", psi_bi);  */
  98.  
  99.     for (knx = 0;  knx < 2;  knx ++)
  100.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  101.             for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  102.                 Psi_M[knx][inx][jnx] = 0.0;
  103. /*    only_on (OBSERVERS) Print_Mat ("Psi_M", Psi_M);  */
  104.   
  105.     for (knx = 0;  knx < 2;  knx ++)
  106.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  107.             for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  108.                 Psi[knx][inx][jnx] = 0.0;
  109. /*    only_on (OBSERVERS) Print_Mat ("Psi", Psi);  */
  110.   
  111.     Compute_Tau ();
  112.     Exchange (Tauz);
  113. /*    only_on (OBSERVERS) Print_Mat ("Tauz", Tauz);  */
  114.  
  115.     for (inx = 0;  inx < num_iter;  inx ++)
  116.         Do_One_Iteration ();
  117.  
  118.     for (inx = 1;  inx < X_PTS;  inx ++)
  119.         for (jnx = 1;  jnx < Y_PTS;  jnx ++)
  120.             Psi_Um[inx][jnx] += Psi_M[0][inx][jnx];
  121. /*    only_on (OBSERVERS) Print_Mat ("Psi_Um", Psi_Um);  */
  122.  
  123.     for (inx = 1;  inx < X_PTS;  inx ++)
  124.         for (jnx = 1;  jnx < Y_PTS;  jnx ++)
  125.             Psi_Lm[inx][jnx] += Psi_M[1][inx][jnx];
  126. /*    only_on (OBSERVERS) Print_Mat ("Psi_Lm", Psi_Lm);  */
  127.  
  128.     only_on (OBSERVER)
  129.         {time = Second_Look (N_TIMER);
  130.          CMMD_send (CMMD_host_node , TIME_TAG, &time, sizeof (double));
  131.         }
  132.    }
  133.  
  134.  
  135. void Compute_My_Coordinates ()
  136.  
  137.    { int my_id;
  138.  
  139. #if 0
  140.      TOT_PROCS =  numtasks();
  141.      X_PROCS =  NPX(TOT_PROCS);
  142.      Y_PROCS = TOT_PROCS/X_PROCS;
  143.      X_PTS   = (X_TOT_PTS / X_PROCS);
  144.      Y_PTS   = (Y_TOT_PTS / Y_PROCS);
  145.      X_PTS_1 =  (X_PTS + 1);
  146.      Y_PTS_1 = (Y_PTS + 1);
  147.      X_PTS_2 = (X_PTS + 2);
  148.      Y_PTS_2 = (Y_PTS + 2);
  149.      ELEMS_Y2= (DB * Y_PTS_2);
  150. #endif
  151.  
  152.      my_id = CMMD_self_address ();
  153.      my_node = my_id - 1;
  154.  
  155.     my_x    = my_node % X_PROCS;
  156.     my_y    = my_node / X_PROCS;
  157.  
  158.     my_west = my_x == 0 ?           UNDEFINED : my_id - 1;
  159.     my_east = my_x == X_PROCS - 1 ? UNDEFINED : my_id + 1;
  160.     my_sout = my_y == 0 ?           UNDEFINED : my_id - X_PROCS;
  161.     my_nort = my_y == Y_PROCS - 1 ? UNDEFINED : my_id + X_PROCS;
  162.  
  163.     lower_x = my_x == 0           ? 2     : 1;
  164.     upper_x = my_x == X_PROCS - 1 ? X_PTS : X_PTS_1;
  165.     lower_y = my_y == 0           ? 2     : 1;
  166.     upper_y = my_y == Y_PROCS - 1 ? Y_PTS : Y_PTS_1;
  167.    }
  168.  
  169.  
  170. void Compute_Constants ()
  171.  
  172.    {pi = 4.0 * atan (1.0);
  173.     t0 = 0.5e-04;
  174.     h_inv  = 1.0 / h;
  175.     h1_inv = 1.0 / h1;
  176.     hh1    = h1 / h;
  177.     hh3    = h3 / h;
  178.     timst  = 2.0 * dtau;
  179.     psi_bi = 0.0;
  180.     ysca   = (double) (X_TOT_PTS - 1) * res;
  181.     ysca_1 = 0.5 * ysca;
  182.     factor = - t0 * pi / ysca_1;
  183.     eig_2  = - h * f0 * f0 / (h1 * h3 * gpr);
  184.     del_sq_inv  = res * res;
  185.     del_x_sq    =   1.0 / (res * res);
  186.     del_y_sq    =   1.0 / (res * res);
  187.     fact_jacob  = - 1.0 / (12.0 * res * res);
  188.     fact_laplac =   1.0 / (res * res);
  189.    }
  190.  
  191.  
  192. void Compute_F ()
  193.  
  194.    {int inx;
  195.  
  196.     for (inx = lower_x - 1;  inx < upper_x + 1;  inx ++)
  197.         F[inx] = f0 + beta * ((double) (my_x * X_PTS + inx - 1) * 
  198.                  res - ysca_1);
  199.    }
  200.  
  201.  
  202. void Initial_Matrix (Matrix, corner_val, border_val, center_val)
  203. double Matrix[][Y_PTS_2], corner_val, border_val, center_val;
  204.  
  205.    {int inx, jnx;
  206.  
  207.     if (my_x == 0)
  208.        {if (my_y == 0)  Matrix[1][1] = corner_val;
  209.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  210.             Matrix[1][jnx] = border_val;
  211.        }
  212.  
  213.     if (my_x == X_PROCS - 1)
  214.        {if (my_y == Y_PROCS - 1)  Matrix[X_PTS][Y_PTS] = corner_val;
  215.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  216.             Matrix[X_PTS][jnx] = border_val;
  217.        }
  218.  
  219.     if (my_y == 0)
  220.        {if (my_x == X_PROCS - 1) Matrix[X_PTS][1] = corner_val;
  221.         for (inx = lower_x;  inx < upper_x;  inx ++)
  222.             Matrix[inx][1] = border_val; 
  223.        }
  224.  
  225.     if (my_y == Y_PROCS - 1)
  226.        {if (my_x == 0)  Matrix[1][Y_PTS] = corner_val;
  227.         for (inx = lower_x;  inx < upper_x;  inx ++)
  228.             Matrix[inx][Y_PTS] = border_val;
  229.        }
  230.  
  231.     for (inx = lower_x;  inx < upper_x;  inx ++)
  232.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  233.             Matrix[inx][jnx] = center_val;
  234.    }
  235.  
  236.  
  237. double Compute_Sum (Matrix, corner_coef, border_coef, center_coef)
  238. double Matrix[][Y_PTS_2], corner_coef, border_coef, center_coef;
  239.  
  240.    {double sum;
  241.     int inx, jnx;
  242.  
  243.     sum = 0.0;
  244.     if (my_x == 0)
  245.        {if (my_y == 0)  sum += corner_coef * Matrix[1][1];
  246.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  247.             sum += border_coef * Matrix[1][jnx];
  248.        }
  249. /*    only_on (OBSERVERS) printf ("sum = %e\n", sum);  */
  250.  
  251.     if (my_x == X_PROCS - 1)
  252.        {if (my_y == Y_PROCS - 1) sum += corner_coef * Matrix[X_PTS][Y_PTS];
  253.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  254.             sum += border_coef * Matrix[X_PTS][jnx];
  255.        }
  256. /*    only_on (OBSERVERS) printf ("sum = %e\n", sum);  */
  257.  
  258.     if (my_y == 0)
  259.        {if (my_x == X_PROCS - 1)  sum += corner_coef * Matrix[X_PTS][1];
  260.         for (inx = lower_x;  inx < upper_x;  inx ++)
  261.             sum += border_coef * Matrix[inx][1];
  262.        }
  263. /*    only_on (OBSERVERS) printf ("sum = %e\n", sum);  */
  264.  
  265.     if (my_y == Y_PROCS - 1)
  266.        {if (my_x == 0)  sum += corner_coef * Matrix[1][Y_PTS];
  267.         for (inx = lower_x;  inx < upper_x;  inx ++)
  268.             sum += border_coef * Matrix[inx][Y_PTS];
  269.        }
  270. /*    only_on (OBSERVERS) printf ("sum = %e\n", sum);  */
  271.  
  272.     for (inx = lower_x;  inx < upper_x;  inx ++)
  273.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  274.             sum += center_coef * Matrix[inx][jnx];
  275. /*    only_on (OBSERVERS) printf ("sum = %e\n", sum);  */
  276.  
  277.     return (sum);
  278.    }
  279.  
  280.  
  281. void Compute_Tau ()
  282.  
  283.    {double curl, tmpd;
  284.     int inx, jnx;
  285.  
  286.     for (inx = lower_x - 1;  inx < upper_x + 1;  inx ++)
  287.       
  288.        {tmpd = pi * (double) (my_x * X_PTS + inx - 1) * 
  289.            res / ysca_1;
  290.     curl = factor * sin (0.0);
  291.         for (jnx = lower_y - 1;  jnx < upper_y + 1;  jnx ++)
  292.             Tauz[inx][jnx] = curl;
  293.        }
  294.    }
  295.  
  296.  
  297. void Do_One_Iteration ()
  298.  
  299.    {int inx, jnx, knx;
  300.     double f4, psi_ai, partial_sum, Compute_Sum ();
  301.     void Laplacian (), Jacobian (), Exchange (), Solve_Helmholtz ();
  302.     void Compute_Ga_and_Gb (), Print_Mat ();
  303.  
  304.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  305.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  306.             Ga[inx][jnx] = 0.0;
  307. /*    only_on (OBSERVERS) Print_Mat ("Ga", Ga);  */
  308.  
  309.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  310.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  311.             Gb[inx][jnx] = 0.0;
  312. /*    only_on (OBSERVERS) Print_Mat ("Gb", Gb);  */
  313.  
  314.     for (knx = 0;  knx < 2;  knx ++)
  315.        {Laplacian (Psi[knx], Work_1[knx]);
  316.     Exchange  (Work_1[knx]);
  317. /*        only_on (OBSERVERS) Print_Mat ("Work_1[knx]", Work_1[knx]);  */
  318.        }
  319.  
  320.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  321.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  322.             Work_2[inx][jnx] = Psi[0][inx][jnx] - Psi[1][inx][jnx];
  323. /*    only_on (OBSERVERS) Print_Mat ("Work_2", Work_2);  */
  324.  
  325.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  326.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  327.             Work_3[inx][jnx] = hh3 * Psi[0][inx][jnx] + hh1 * Psi[1][inx][jnx];
  328. /*    only_on (OBSERVERS) Print_Mat ("Work_3", Work_3);  */
  329.  
  330.     for (knx = 0;  knx < 2;  knx ++)
  331.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  332.             for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  333.                 Temp[knx][inx][jnx] = Psi[knx][inx][jnx];
  334. /*    only_on (OBSERVERS) Print_Mat ("Temp", Temp);  */
  335.  
  336.     for (knx = 0;  knx < 2;  knx ++)
  337.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  338.             for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  339.                 Psi[knx][inx][jnx] = Psi_M[knx][inx][jnx];
  340. /*    only_on (OBSERVERS) Print_Mat ("Psi[0]", Psi[0]);  */
  341. /*    only_on (OBSERVERS) Print_Mat ("Psi[1]", Psi[1]);  */
  342.  
  343.     for (knx = 0;  knx < 2;  knx ++)
  344.        {Laplacian (Psi_M[knx], Work_7[knx]);
  345.     Exchange  (Work_7[knx]);
  346. /*        only_on (OBSERVERS) Print_Mat ("Work_7[knx]", Work_7[knx]);  */
  347.        }
  348.  
  349.     for (knx = 0;  knx < 2;  knx ++)
  350.         for (inx = lower_x - 1;  inx < upper_x + 1;  inx ++)
  351.             for (jnx = lower_y - 1;  jnx < upper_y + 1;  jnx ++)
  352.                 Work_1[knx][inx][jnx] += F[inx];
  353. /*    only_on (OBSERVERS) Print_Mat ("Work_1[0]", Work_1[0]);  */
  354. /*    only_on (OBSERVERS) Print_Mat ("Work_1[1]", Work_1[1]);  */
  355.  
  356.     for (knx = 0;  knx < 2;  knx ++)
  357.        {Jacobian (Work_1[knx], Temp[knx], Work_5[knx]);
  358.     Exchange (Work_5[knx]);
  359. /*        only_on (OBSERVERS) Print_Mat ("Work_5[knx]", Work_5[knx]);  */
  360.        }
  361.  
  362.     for (knx = 0;  knx < 2;  knx ++)
  363.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  364.             for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  365.                 Psi_M[knx][inx][jnx] = Temp[knx][inx][jnx];
  366. /*    only_on (OBSERVERS) Print_Mat ("Psi_M[0]", Psi_M[0]);  */
  367. /*    only_on (OBSERVERS) Print_Mat ("Psi_M[1]", Psi_M[1]);  */
  368.  
  369.     for (knx = 0;  knx < 2;  knx ++)
  370.        {Laplacian (Work_7[knx], Work_4[knx]);
  371.     Exchange  (Work_4[knx]);
  372. /*        only_on (OBSERVERS) Print_Mat ("Work_4[knx]", Work_4[knx]);  */
  373.        }
  374.  
  375.     Jacobian (Work_2, Work_3, Work_6);
  376.     Exchange (Work_6);
  377. /*    only_on (OBSERVERS) Print_Mat ("Work_6", Work_6);  */
  378.  
  379.     for (knx = 0;  knx < 2;  knx ++)
  380.        {Laplacian (Work_4[knx], Work_7[knx]);
  381.     Exchange  (Work_7[knx]);
  382. /*        only_on (OBSERVERS) Print_Mat ("Work_7[knx]", Work_7[knx]);  */
  383.        }
  384.  
  385.     Compute_Ga_and_Gb ();
  386.     Exchange (Ga);
  387.     Exchange (Gb);
  388. /*    only_on (OBSERVERS) Print_Mat ("Ga", Ga);  */
  389. /*    only_on (OBSERVERS) Print_Mat ("Gb", Gb);  */
  390.  
  391.     Solve_Helmholtz (eig_2, Ga, Old_Ga);
  392. /*    only_on (OBSERVERS) Print_Mat ("Ga", Ga);  */
  393.  
  394.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  395.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  396.             Old_Ga[inx][jnx] = Ga[inx][jnx];
  397. /*    only_on (OBSERVERS) Print_Mat ("Old_Ga", Old_Ga);  */
  398.  
  399.     partial_sum = Compute_Sum (Ga, 0.25, 0.50, 1.00);
  400.     psi_ai = CMMD_reduce_double (partial_sum, CMMD_combiner_add);
  401. /*    only_on (OBSERVERS) printf ("psi_ai = %e\n", psi_ai);  */
  402.  
  403.     f4 = - psi_ai / psi_bi;
  404.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  405.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  406.             Ga[inx][jnx] += f4 * Psi_B[inx][jnx];
  407. /*    only_on (OBSERVERS) Print_Mat ("Ga", Ga);  */
  408.  
  409.     Solve_Helmholtz (eig_2, Gb, Old_Gb);
  410. /*    only_on (OBSERVERS) Print_Mat ("Gb", Gb);  */
  411.  
  412.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  413.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  414.             Old_Gb[inx][jnx] = Gb[inx][jnx];
  415. /*    only_on (OBSERVERS) Print_Mat ("Old_Gb", Old_Gb);  */
  416.  
  417.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  418.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  419.             {Work_2[inx][jnx] = Gb[inx][jnx] - hh1 * Ga[inx][jnx];
  420.              Work_3[inx][jnx] = Gb[inx][jnx] + hh3 * Ga[inx][jnx];
  421.         }
  422. /*    only_on (OBSERVERS) Print_Mat ("Work_2", Work_2);  */
  423. /*    only_on (OBSERVERS) Print_Mat ("Work_3", Work_3);  */
  424.  
  425.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  426.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  427.             Psi[0][inx][jnx] += timst * Work_3[inx][jnx];
  428. /*    only_on (OBSERVERS) Print_Mat ("Psi[0]", Psi[0]);  */
  429.  
  430.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  431.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  432.             Psi[1][inx][jnx] += timst * Work_2[inx][jnx];
  433. /*    only_on (OBSERVERS) Print_Mat ("Psi[1]", Psi[1]);  */
  434.    }
  435.  
  436.  
  437. void Exchange (Mat)
  438. double Mat[][Y_PTS_2];
  439.  
  440.    {if (X_PROCS > 1)
  441.         if (my_x != 0 && my_x != X_PROCS - 1)
  442.            {Send_Receive (my_west, &Mat[0][0], my_east, &Mat[X_PTS][0]);
  443.             Send_Receive (my_east, &Mat[X_PTS_1][0], my_west, &Mat[1][0]);
  444.            }
  445.         else
  446.            if (my_x == 0)
  447.               {Send    (my_east, &Mat[X_PTS][0]);
  448.                Receive (my_east, &Mat[X_PTS_1][0]);
  449.           }
  450.            else
  451.               {Receive (my_west, &Mat[0][0]);
  452.                Send    (my_west, &Mat[1][0]);
  453.           }
  454.  
  455.     if (Y_PROCS > 1)
  456.         if (my_y != 0 && my_y != Y_PROCS - 1)
  457.            {Send_Receive_V (my_sout, &Mat[0][0], my_nort, &Mat[0][Y_PTS]);
  458.             Send_Receive_V (my_nort, &Mat[0][Y_PTS_1], my_sout, &Mat[0][1]);
  459.            }
  460.         else
  461.            if (my_y == 0)
  462.               {Send_V    (my_nort, &Mat[0][Y_PTS]);
  463.                Receive_V (my_nort, &Mat[0][Y_PTS_1]);
  464.           }
  465.            else
  466.               {Receive_V (my_sout, &Mat[0][0]);
  467.                Send_V    (my_sout, &Mat[0][1]);
  468.           }
  469.    }
  470.  
  471.  
  472. void Compute_Ga_and_Gb ()
  473.  
  474.    {int inx, jnx;
  475.  
  476.     if (my_x == 0)
  477.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++)
  478.            {Ga[1][jnx] = Work_5[0][1][jnx] + lf * Work_7[0][1][jnx] - lf *
  479.                          Work_7[1][1][jnx];
  480.  
  481.             Gb[1][jnx] = hh1 * Work_5[0][1][jnx] + hh3 * 
  482.                                Work_5[1][1][jnx] + lf * hh1 * 
  483.                                Work_7[0][1][jnx] + lf * hh3 *
  484.                                Work_7[1][1][jnx];
  485.        }
  486.  
  487.     if (my_x == X_PROCS - 1)
  488.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++)
  489.            {Ga[X_PTS][jnx] = Work_5[0][X_PTS][jnx] + lf * 
  490.                              Work_7[0][X_PTS][jnx] - lf *
  491.                              Work_7[1][X_PTS][jnx];
  492.  
  493.             Gb[X_PTS][jnx] = hh1 * Work_5[0][X_PTS][jnx] + hh3 * 
  494.                                    Work_5[1][X_PTS][jnx] + lf * hh1 * 
  495.                                    Work_7[0][X_PTS][jnx] + lf * hh3 *
  496.                                    Work_7[1][X_PTS][jnx];
  497.        }
  498.  
  499.     if (my_y == 0)
  500.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  501.            {Ga[inx][1] = Work_5[0][inx][1] + lf * Work_7[0][inx][1] - lf * 
  502.                          Work_7[1][inx][1];
  503.  
  504.             Gb[inx][1] = hh1 * Work_5[0][inx][1] + hh3 * 
  505.                                Work_5[1][inx][1] + lf * hh1 * 
  506.                                Work_7[0][inx][1] + lf * hh3 * 
  507.                                Work_7[1][inx][1];
  508.        }
  509.  
  510.     if (my_y == Y_PROCS - 1)
  511.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  512.            {Ga[inx][Y_PTS] = Work_5[0][inx][Y_PTS] + lf * 
  513.                              Work_7[0][inx][Y_PTS] - lf * 
  514.                              Work_7[1][inx][Y_PTS];
  515.  
  516.             Gb[inx][Y_PTS] = hh1 * Work_5[0][inx][Y_PTS] + hh3 * 
  517.                                    Work_5[1][inx][Y_PTS] + lf * hh1 * 
  518.                                    Work_7[0][inx][Y_PTS] + lf * hh3 * 
  519.                                    Work_7[1][inx][Y_PTS];
  520.        }
  521.  
  522.     for (inx = lower_x;  inx < upper_x;  inx ++)
  523.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  524.            {Ga[inx][jnx] = Work_5[0][inx][jnx] - Work_5[1][inx][jnx] + eig_2 *
  525.                            Work_6[inx][jnx] + h1_inv * Tauz[inx][jnx] + lf *
  526.                            Work_7[0][inx][jnx] - lf * Work_7[1][inx][jnx];
  527.  
  528.             Gb[inx][jnx] = hh1 * Work_5[0][inx][jnx] + hh3 * 
  529.                            Work_5[1][inx][jnx] + h_inv * Tauz[inx][jnx] + lf *
  530.                            hh1 * Work_7[0][inx][jnx] + lf * hh3 *
  531.                            Work_7[1][inx][jnx];
  532.        }
  533.    }
  534.  
  535.  
  536. void Solve_Helmholtz (elmbda, Mat_Out, Mat_Ini)
  537. double Mat_Out[][Y_PTS_2], Mat_Ini[][Y_PTS_2];
  538. double elmbda;
  539.  
  540.    {int inx, jnx, knx, lower_ch_y, lower_top, offset, my_flag, all_flag;
  541.     double old_val, fac;
  542.     void Print_Mat ();
  543.  
  544.     fac = 1.0 / (4.0 - del_sq_inv * elmbda);
  545.     if (my_x == 0)
  546.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++)
  547.             Work_8[1][jnx] = Mat_Ini[1][jnx];
  548.  
  549.     if (my_x == X_PROCS - 1)
  550.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++)
  551.             Work_8[X_PTS][jnx] = Mat_Ini[X_PTS][jnx];
  552.  
  553.     if (my_y == 0)
  554.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  555.             Work_8[inx][1] = Mat_Ini[inx][1];
  556.  
  557.     if (my_y == Y_PROCS - 1)
  558.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  559.             Work_8[inx][Y_PTS] = Mat_Ini[inx][Y_PTS];
  560.  
  561.     for (inx = lower_x;  inx < upper_x;  inx ++)
  562.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  563.             Work_8[inx][jnx] = omega * fac * (Mat_Ini[inx][jnx+1] +
  564.                         Mat_Ini[inx][jnx-1] + Mat_Ini[inx+1][jnx] + 
  565.                         Mat_Ini[inx-1][jnx] - del_sq_inv * Mat_Out[inx][jnx]) +
  566.                         (1 - omega) * Mat_Ini[inx][jnx];
  567.     Exchange (Work_8);
  568. /*    only_on (OBSERVERS) Print_Mat ("Work_8", Work_8);  */
  569.  
  570.     do {my_flag = FALSE;
  571.     for (knx = 0;  knx < 2;  knx ++)
  572.       {lower_ch_y = lower_y +  (lower_x + lower_y + knx) % 2;
  573.            lower_top  = 1 -  2  * ((lower_x + lower_y + knx) % 2);
  574.             for (inx = lower_x, offset = 0;  inx < upper_x;  
  575.                 inx ++, offset = lower_top - offset)
  576.                for (jnx = lower_ch_y + offset;  jnx < upper_y;  jnx += 2)
  577.                    {old_val = Work_8[inx][jnx];
  578.                     Work_8[inx][jnx] = omega * fac * (Work_8[inx][jnx+1] +
  579.                          Work_8[inx][jnx-1] + Work_8[inx+1][jnx] + 
  580.                          Work_8[inx-1][jnx] - del_sq_inv * Mat_Out[inx][jnx]) +
  581.                          (1 - omega) * Work_8[inx][jnx];
  582.                     if ((old_val - Work_8[inx][jnx] >   tolerance) ||
  583.                         (old_val - Work_8[inx][jnx] < - tolerance)) 
  584.                         my_flag = TRUE;
  585.                 }
  586.             Exchange (Work_8);
  587.            }
  588.  
  589.         all_flag = CMMD_reduce_int (my_flag, CMMD_combiner_ior);
  590. /*        only_on (OBSERVERS) Print_Mat ("Work_8", Work_8);  */
  591.        } while (all_flag == TRUE);
  592.  
  593.     for (inx = 0;  inx < X_PTS_2;  inx ++)
  594.         for (jnx = 0;  jnx < Y_PTS_2;  jnx ++)
  595.             Mat_Out[inx][jnx] = Work_8[inx][jnx];
  596. /*    only_on (OBSERVERS) Print_Mat ("Mat_Out", Mat_Out);  */
  597.    }
  598.  
  599.  
  600. void Laplacian (In, Out)
  601. double In[][Y_PTS_2], Out[][Y_PTS_2];
  602.  
  603.     {int inx, jnx;
  604.  
  605.     if (my_x == 0)       
  606.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++) 
  607.             Out[1][jnx]     = 0.0;
  608.  
  609.     if (my_x == X_PROCS - 1)       
  610.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++) 
  611.             Out[X_PTS][jnx] = 0.0;
  612.  
  613.     if (my_y == 0)       
  614.         for (inx = 1;  inx < X_PTS_1;  inx ++) 
  615.             Out[inx][1]     = 0.0;
  616.  
  617.     if (my_y == Y_PROCS - 1)       
  618.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  619.             Out[inx][Y_PTS] = 0.0;
  620.  
  621.     for (inx = lower_x;  inx < upper_x;  inx ++)
  622.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  623.             Out[inx][jnx] = fact_laplac * (In[inx][jnx+1] + In[inx][jnx-1] + 
  624.                             In[inx+1][jnx] + In[inx-1][jnx] - 4.0 * 
  625.                             In[inx][jnx]);
  626.     }
  627.  
  628.  
  629. void Jacobian (In_1, In_2, Out)
  630. double In_1[][Y_PTS_2], In_2[][Y_PTS_2], Out[][Y_PTS_2];
  631.  
  632.    {double f_1, f_2, f_3, f_4, f_5, f_6, f_7, f_8;
  633.     int inx, jnx;
  634.  
  635.     if (my_x == 0)       
  636.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++) 
  637.             Out[1][jnx]     = 0.0;
  638.  
  639.     if (my_x == X_PROCS - 1)       
  640.         for (jnx = 1;  jnx < Y_PTS_1;  jnx ++) 
  641.             Out[X_PTS][jnx] = 0.0;
  642.  
  643.     if (my_y == 0)       
  644.         for (inx = 1;  inx < X_PTS_1;  inx ++) 
  645.             Out[inx][1]     = 0.0;
  646.  
  647.     if (my_y == Y_PROCS - 1)       
  648.         for (inx = 1;  inx < X_PTS_1;  inx ++)
  649.             Out[inx][Y_PTS] = 0.0;
  650.  
  651.     for (inx = lower_x;  inx < upper_x;  inx ++)
  652.         for (jnx = lower_y;  jnx < upper_y;  jnx ++)
  653.            {f_1 = (In_2[inx-1][jnx] + In_2[inx-1][jnx+1] - In_2[inx+1][jnx] -
  654.                    In_2[inx+1][jnx+1]) * (In_1[inx][jnx+1] - In_1[inx][jnx]);
  655.  
  656.             f_2 = (In_2[inx-1][jnx-1] + In_2[inx-1][jnx] - In_2[inx+1][jnx-1] -
  657.                    In_2[inx+1][jnx]) * (In_1[inx][jnx] - In_1[inx][jnx-1]);
  658.  
  659.             f_3 = (In_2[inx][jnx+1] + In_2[inx+1][jnx+1] - In_2[inx][jnx-1] -
  660.                    In_2[inx+1][jnx-1]) * (In_1[inx+1][jnx] - In_1[inx][jnx]);
  661.  
  662.             f_4 = (In_2[inx-1][jnx+1] + In_2[inx][jnx+1] - In_2[inx-1][jnx-1] -
  663.                    In_2[inx][jnx-1]) * (In_1[inx][jnx] - In_1[inx-1][jnx]);
  664.  
  665.             f_5 = (In_2[inx][jnx+1]   - In_2[inx+1][jnx]) * 
  666.                   (In_1[inx+1][jnx+1] - In_1[inx][jnx]);
  667.  
  668.             f_6 = (In_2[inx-1][jnx]   - In_2[inx][jnx-1]) * 
  669.                   (In_1[inx][jnx]     - In_1[inx-1][jnx-1]);
  670.  
  671.             f_7 = (In_2[inx+1][jnx]   - In_2[inx][jnx-1]) * 
  672.                   (In_1[inx+1][jnx-1] - In_1[inx][jnx]);
  673.  
  674.             f_8 = (In_2[inx][jnx+1]   - In_2[inx-1][jnx]) * 
  675.                   (In_1[inx][jnx]     - In_1[inx-1][jnx+1]);
  676.  
  677.             Out[inx][jnx] = fact_jacob * (f_1 + f_2 + f_3 + f_4 + 
  678.                                           f_5 + f_6 + f_7 + f_8);
  679.            }
  680.     }
  681.  
  682.  
  683. void Print_Mat (name, Matrix)
  684. char *name;
  685. double Matrix[][Y_PTS_2];
  686.  
  687.    {int inx, jnx, dummy;
  688.  
  689.     if (TOT_PROCS > 1 && my_node != 0) 
  690.         CMMD_receive (my_node - 1, ANY_TAG, &dummy, sizeof (int));
  691.  
  692.     printf ("\n");
  693.     if (my_node == 0) printf ("            %s\n", name);
  694.     if (my_node == 0) printf ("            --------\n");
  695.     for (jnx = Y_PTS_1;  jnx >= 0;  jnx --)
  696.        {printf ("[%02d]", jnx);
  697.         for (inx = 0;  inx < X_PTS_2;  inx ++)
  698.             printf (" %9.2e", Matrix[inx][jnx]);
  699.         printf ("\n");
  700.        }
  701.  
  702.     if (TOT_PROCS > 1)
  703.         if (my_node != 0) 
  704.             CMMD_send ((my_node+1)%TOT_PROCS, ANY_TAG, &dummy, sizeof (int));
  705.         else
  706.            {CMMD_send    (1,             ANY_TAG, &dummy, sizeof (int));
  707.             CMMD_receive (TOT_PROCS - 1, ANY_TAG, &dummy, sizeof (int));
  708.            }
  709.    }
  710.